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1 .  INTRODUCTION 


Various  mathematical  models  exist  for  the  purpose  of  re¬ 
presenting  the  anomalous  gravity  field  of  the  Earth:  integral 
formulas,  series  of  harmonic  functions,  etc.  Due  to  the  harmonicitv 
of  the  underlying  anomalous  potential  outside  the  Earth's  surface, 
the  integral  kernels  and  the  solid  spherical  harmonics  are  different 
from  zero  almost  everywhere;  spherical  harmonic  representations 
have,  in  addition,  the  following  disadvantages:  the  solution  reacts 
globally  to  local  changes;  in  order  to  represent  the  short- 
periodic  part  of  the  spectrum,  the  series  development  would  require 
a  very  high  degree;  although  spherical  harmonics  are  evaluated  by 
the  use  of  recursion  formulas,  a  series  development  of  that  kind 
remains  an  expensive  task.  Low  degree  spherical  harmonic  coefficients 
are  known  to  permit  a  simple  physical  interpretation  in  terms  of 
mass,  coordinates  of  the  center  of  gravity,  moments  of  inertia, 
etc.;  for  higher  degrees  no  simple  physical  interpretation  is 
known.  Its  determination  is  possible  without  knowing  anything 
about  the  mass  distribution  inside  the  Earth.  The  price  that  has 
to  be  paid  for  this  ignorance  is  the  global  behavior.  Neither 
the  numerical  evaluation  of  integral  formulas  nor  the  evaluation 
of  a  long  harmonic  series  are  adequate  means  for  fast  (real-time) 
and  flexible  (local)  representations  of  the  gravity  field.  The 
cause  of  the  gravity  field,  the  mass  distribution  inside  the  Earth, 
is  not  transparent  either  in  such  kind  of  representations. 

Recently  two  virtually  quite  different  gravity  field  repre¬ 
sentation  techniques  are  being  developed,  which  will  possibly 
replace  (but  for  sure  supplement)  existing  techniques: 
a)  The  finite  element  representation  of  the  gravity  field  outside 
the  Earth.  This  method  is  extremely  powerful  as  far  as  the  fast 
prediction  of  gravity  field  quantities  is  concerned;  it  is  based 
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on  functions  with  local  support  and  therefore,  the  represented 
anomalous  potential  is  not  strictly  harmonic;  a  combination  with 
spherical  harmonic  solutions  seems  to  be  very  difficult;  there 
is  practically  no  relation  to  the  mass  distribution  inside  the 
Earth. 

b)  The  mass  model  representation  of  the  anomalous  gravity  field 
which,  at  the  present  state  of  the  art,  is  restricted  to  its 
simplest  form,  the  point  mass  models.  This  approach  towards  the 
cause  of  the  gravity  field  is  at  the  same  time  an  approach  of 
geodesy  and  geophysics.  The  relation  of  the  model  to  the  physical 
reality  is  transparent,  geophysical  evidencies  in  terms  of  geo¬ 
logical  structures  can  be  easily  implemented  into  the  model.  The 
harmonicity  of  the  represented  anomalous  potential  is  guaranteed. 
The  reciprocal  distance  and  its  derivatives  provide  the  simplest 
possible  relation  between  the  data  (point  masses)  and  any  derived 
gravity  field  quantity.  Due  to  the  global  support  of  the  reci¬ 
procal  distance,  all  point  masses  contribute  to  the  prediction  of 
a  single  quantity;  however,  the  remote  zone  effect  can  be  used 
advantageously  to  reduce  the  number  of  actually  used  point  masses. 
The  relation  between  point  masses  and  harmonic  coefficients  of 
the  anomalous  gravitational  potential  is  straightforward  and 
simple . 

There  is  of  course  an  essential  drawback  in  the  application 
of  this  technique:  theoretically  there  exists  an  infinite  number 
of  mass  distributions  which  are  compatible  with  the  observed 
gravity  field.  From  geophysical  evidences  we  know  the  main  features 
of  the  density  distribution.  We  can  make  a  virtue  of  necessity 
and  select  a  solution  which  is  both  simple  and  geophysically  re¬ 
levant.  As  far  as  the  distribution  of  point  masses  is  concerned, 
we  strongly  emphasize  the  importance  of  regular  patterns.  A  proper 
design  of  the  data  pattern  can  reduce  the  calculation  efforts 
dramatically  if  the  algorithm  is  sophisticated  enough  to  realize 
the  data  geometry. 
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Point  mass  models  wi^l  be  determined  primarily  on  the  basis 
of  harmonic  coefficients  and/or  surface  mean  gravity  anomalies. 

In  order  to  get  an  idea  about  the  relation  between  these  quantities, 
we  have  to  start  at  the  very  beginning  and  investigate  the  re¬ 
sponse  of  essential  statistical  quantities  of  the  anomalous  gravi¬ 
tational  field  to  changes  in  the  spatial  point  mass  model  arrange¬ 
ment.  Statistical  models  of  continuous  mass  distributions  provide 
important  insight  from  a  theoretical  point  of  view,  an  investigation 
of  models  of  discrete  distributions  in  terms  of  point  masses  is 
indispensable  from  a  practical  point  of  view.  In  the  latter  case 
the  enormously  powerful  tool  of  frequency  domain  methods  on  the 
sphere  is  a  prerequisite  for  serious  model  calculations.  The 
variation  of  the  variance  and  correlation  length  of  a  homogeneous 
and  isotropic  gravity  anomaly  covariance  function,  generated  by 
mass  distributions  in  various  depths,  is  discussed  in  detail.  For 
discrete  distributions  the  relation  to  harmonic  coefficients  of 
the  anomalous  gravitational  potential  is  studied;  guidelines  for 
actual  computations,  using  the  concept  of  Fast  Fourier  Transform 
on  the  sphere,  are  given.  Since  a  mass  model  should  reproduce  as 
close  as  possible  the  power  per  degree  of  the  anomalous  gravitational 
field,  emphasis  is  put  on  the  numerical  estimation  of  degree 
variances.  The  actual  calculation  of  a  point  mass  model  from  real 
world  data  will  be  the  subject  of  a  forthcoming  report. 
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2.  DEPTHS  OF  ANOMALOUS  MASS  DISTRIBUTIONS  AND  THE  GRAVITY 
ANOMALY  COVARIANCE  FUNCTION 


The  essential  statistical  properties  of  the  gravity 
anomaly  field  can  be  described  by  a  few  parameters,  the  variance 
C0  ,  the  correlation  length  5  ,  and  the  variance  of  the  horizontal 
gravity  gradient  ,  Go  .  The  latter  shows  a  strong  response  to 
topographical  and  local  density  anomalies;  the  correlation  length 
and  the  variance  are  much  less  affected  by  local  anomalies  and 
should  therefore  respond  to  deeper  density  anomalies.  It  is  of 
primary  concern  to  know  the  dependence  of  these  two  quantities 
on  the  depth  of  the  mass  point  level  D  . 

In  order  to  investigate  this  problem,  let  us  start  with  a 
highly  unlikely  but  nevertheless  very  instructive  case,  we  assume 
a  continuous  mass  distribution  at  the  depth  D  below  the  surface 
of  the  mean  terrestrial  sphere,  with  zero  average  (positive  and 
negative  density  anomalies) ,  and  a  white  noise  covariance  function; 
the  center  of  gravity  is  supposed  to  coincide  with  the  origin  of 
the  coordinate  system.  Needless  to  say,  this  model  is  far  from 
reality;  it  is  a  very  pessimistic  model  insofar  as  the  anomalous 
masses  are  assumed  to  be  uncorrelated.  In  the  following  we  shall 
investigate,  how  the  gravity  anomaly  covariance  function  at  mean 
sea  level  responds  to  that  white  noise  mass  distribution  at  depth 
D. 

The  disturbing  potential  T  is  given  by  (Heiskanen  and 
Moritz,  1967,  p.  5) 
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T(p)  -  G  j'jr  tt§W  do<o> 


with  G  .  gravitational  constant, 

l  .  distance  (P,Q) 

do  . ■  element  of  solid  angle 

udo  .  element  of  anomalous  mass. 


(2.1) 


With 

l2  =  r2  +  (R-D)2  -  2r(R-D)cos* 

and 

Iy  =  r  -  (R-D)  cos  ip  , 

the  radial  derivate  of  T  at  r  =  R  is  obtained  by 


rr  [R-  (R-  D)  cosi))  J  u(0) 

=  -G  j  - ^ -  do  (Q )  (2.2) 

'r  =  R  J  1  e3(  P.Q) 

0 

with 

R  ...  mean  earth  radius, 

D  ...  depth  of  mass  anomaly  layer, 

'■P  ...  spherical  distance 

Both  (2.1)  and  (2.2)  describe  a  linear  system  with  input  u  and 
output  T  and  3T/3r  ,  respectively;  the  integral  kernel  Gi"1 
and  -G[R-  (R-D)  cosip  ]  S.-3  ,  resp.,  is  the  corresponding  system's 


3T(P) 

3r 
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impulse  response.  Due  to  ag  =  -3T/3r  -  2T/r  ,  the  impulse 
response  of  the  linear  system  with  input  u  and  output  ag 
is  equal  to  Gl~i{R-  (R-D)cos y  -  2l2/R)  . 

In  order  to  derive  the  gravity  anomaly  covariance  function 
C  from  the  (white  noise)  mass  anomaly  covariance  function  C 
we  have  to  know  the  system  "function"  which  is  nothing  else  chan 
the  infinite  vector  of  the  integral  kernel's  eigenvalues.  The 
eigenvalues  \  of  the  isotropic  integral  kernel  Gz_5[R- 
-  (R-D)  cost*/  -  2i  VR  i  equal  its  projection  onto  the  set  of  Le¬ 
gendre  polynomials,  multiplied  by  2n  , 

1 

r 

\n  =2-  j  G*”3[R-  (R-  D)t  -  222/R.'Pn  (t)dt  (2.3) 

-  1 

with  t  =  cosy  and  P  denoting  the  Legendre  polynomial  of 
degree  n  (Muller,  1966,  p.  20).  The  integral  can  be  solved 
for  arbitrary  n  by  representing  i~ 1  and  Z~3  in  terms  of  a 
series  of  Legendre  polynomials,  and  observing  the  well-known  re¬ 
currence  and  orthogonality  relations.  A  much  simpler  way  is  to 
consider  the  linear  system  with  input  u  and  output  ig  as  a 
linear  system  in  cascade,  consisting  of  two  linear  subsystems; 
the  output  of  the  first  system  is  input  to  the  second  one.  Here 
the  first  system  is  represented  by  equation  (2.1),  the  second 
system  is  simply  the  boundary  condition 

=  -21  2T  .  (2.4) 

ar  r 

The  cascade  system  function  equals  the  product  of  the  individual 
system  functions  (Papoulis,  1968,  pp.  50,  51).  System  1  (eq.  2.1) 
has  eigenvalues 


V  i  i 


(2.5) 


=2 ?_1P  (t)dt  ; 
J  n 
.  1 


expressing  i  _1  in  terms  of  the  series 


n  "  i  —  1  Y  P 

R  -  n 
n  =  0 


with  oi :  =  1  -  -p  ,  and  observing  the  orthogonality  relations  of 


Legendre  polynomials  Pm>  =  6  n  m2/(2n  +  1  ^  '  we  obtain 


(  1 )  4  tG  a 1 


R(2n  +  I! 


.2.5) 


System  2  (eq.  2.4)  has  eigenvalues 


n  -  1 


(2.6) 


(Heiskanen  and  Moritz,  1967,  p.  97);  therefore  the  cascade  eigen 


values  (Cascade  system  "function";  are  given  by  <  =  ( 1 !  .  ■.  2) 

an  n 


4^G 


_  n  -  1  „  n 

\  —  _  ■  •  - »■  Cl 

’a  R2  2n  +  1 


:2.7: 


We  nave  assumed  that  the  mass  anomalies  are  uncorrelated 
(white  noise) ;  the  corresponding  covariance  function  is 
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with  the  variance  M0  and  the  Dirac  distribution  6(t)  .  Expres¬ 
sing  C  ,  in  terms  of  a  series  of  Legendre  polynomials, 


C  (t)  =  i  v  P  (t) 
u  u  -  n  ri 


we  obtain  with  (2-8)  the  degree  variances  y 


M0  :  5 ( t  -  1 ) P  ( t) dt 


2n+l 


which  reduces,  due  to  the  reproducing  property  of  6 ( t)  ,  to 


y  =  ~  ( 2n  +  i; 
n  4 


(2.9) 


(y  }  is  the  power  spectrum  of  the  white  noise  anomalous  mass 
distribution.  The  power  spectrum  {grj}  of  the  gravity  anomaly 
field  is  then  given  by 


gn  = 


gn  = 


f  4ttG]2  Me  ( n- 1 ) 2  2n 

lTP“J  2  2n  +  1  a 


(2.10) 


and  the  gravity  anomaly  covariance  function  by 


C  (t)  ■  l  g  P  (t)  . 

gg  l2  ’n  n 


(2.11) 


The  behavior  of  a2n (n  -  1) 2/ (2n +  1}  is  graphically  represented 
in  Fig.  2.1a,b.  The  area  below  the  curves  is  very  closely  rela¬ 
ted  to  the  total  variance  C  (0)  ;  it  is  quite  obvious  that 
the  variance  drops  dramatically  with  increasing  deptn  D  .  (Note 
that  for  all  curves  a  common  factor  (4ttG/R2)2-  M0/2  =  1  has  been 
used.)  The  maximum  contribution  to  the  variance  decreases  towards 
the  low  degree  with  increasing  depth. 


1 


.1 


1 


I 

1 

j 


s 

tM 


40 


50  km 
30  km 


degree  n 


Fig. 


2.1a,b  Gravity  anomaly  degree  variances  due 

to  a  white  noise  anomalous  mass  distri¬ 
bution  at  various  depths  D  ; 

(4itG/R2)2  •  M0/2  =  1  has  been  assumed. 


The  simple  form  of  the  gravity  anomaly  degree  variances  (2.10) 
allows  the  variance  to  be  represented  by  a  closed  expression. 


1 1 


Observing 


(n-1 ) 2_  n  59  1 

2n  +  1  ~  2  4  4  2n  +  1  ' 


*  i 

V  „2n  _  X _ 

l  a 


i  ~  i 


a  <  1 


n*0 


(Ryshik  and  Gradstein,  1963,  No.  1.231,  p.  24)  and  noting  that 


n=0 


r,„2n  _  a.  d 
n“  ‘  2  da 


r  2  n 


n  =  0 


( 1  “  V  )' 


i03^r  =  ;  j  -  £  ln  Hi  • 

n=0  ‘  n=0 


we  obtain  the  closed  expression 


Ca  (0)  ■=  M0(if)2 


9  ,  1  +  a  .  14a2  -  10 

—  in  - -  +  ~r\ - rTT  0 

a  1-a  ( 1  -  a  “  )  *• 


(2.12] 


For  moderate  values  of  D  (DS  100  km)  ,  the  variance  C  (0) 

y  y 

easily  be  shown  to  be  approximately  proportional  to  D'2  , 


c„to>  ■  ^ 


(2.i3: 


can 


With  other  words,  the  gravity  anomaly  variance  due  to  a  white 
noise  anomalous  mass  distribution  at  depth  D  decreases  with  the 


COVflRIfiNCE  (DG.DG),  NORMAL 1  ZED 
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square  of  the  depth  -  a  very  remarkable  result.  Consequently,  approx¬ 
imately  the  same  gravity  anomaly  variance  at  zero  level  is  genera¬ 
ted  by  two  anomalous  mass  distributions  at  level  Dlf  and  Q-  , 
if  the  corresponding  mass  anomaly  variances  Mi  and  M2  benave 
like 


M; 

M: 


Therefore,  we  conclude  that  the  ratios  M,/D2  must  be  selected 

i  i 

properly,  if  the  gravity  anomaly  field's  power  is  to  be  reproduced 
by  mass  anomaly  distributions  at  various  depths . 

Fig.  2.2a,b  show  the  actual  gravity  anomaly  (normalized)  co- 
variance  functions,  produced  by  white  noise  anomalous  mass  distri¬ 
butions.  It  is  obvious  that  the  correlation  length  5  increases 


DISTANCE  1KM) 


COVARIANCE  (DG.DG).  NORMALIZED 
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DISTANCE  (KM) 


Fig.  2.2a,b  Gravity  anomaly  covariance  functions 
at  sea  level  due  to  a  white  noise 
anomalous  mass  distribution  at  various 
depths  D  ;  the  covariances  are  normali¬ 
zed  to  C  (0)  =  1  . 


with  increasing  depth  D  ;  it  is  less  obvious  that  $  depends 
almost  linearly  on  D  (at  least  for  moderate  values  of  D)  with 
a  proportionality  factor  close  to  3/2  , 
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Z  -  |  D  .  (2.14) 

From  these  figures  we  conclude  that  the  global  gravity  anomaly  co- 
variance  function  with  a  correlation  length  of  about  45  km  cannot 
be  generated  by  mass  anomalies,  located  below  30  km  ,  alone;  shal¬ 
low  mass  anomalies  must  considerably  contribute  to  the  observed 
45  km  correlation  length.  The  topographic  masses  and  near  sur¬ 
face  mass  anomalies  alone,  vice  versa,  can  hardly  account  for  a 
correlation  length  of  45  km.  The  problem  of  relating  the  known 
gravity  anomaly  field  to  unknown  mass  anomalies  is  a  (difficult) 
matter  of  tuning  our  sensors  to  proper  frequencies.  This  delicate 
situation  has  been  described  best  by  Alfred  Wegener  (1929): 

"We  are  like  a  judge  confronted  by  a  defendent  who  declines 
to  answer,  and  we  must  determine  the  truth  from  the  circumstancial 
evidence" . 


3.  DISCRETE  POINT  MASSES  AND  CORRESPONDING  COVARIANCE 
FUNCTIONS 

An  anomalous  mass  distribution  represented  by  a  white  noise 
covariance  function  as  investigated  in  the  foregoing  chapter, 
provides  insight  into  the  relation  between  anomalous  mass  and 
anomalous  gravity.  That  highly  artificial  case  is  quite  unlikely 
to  be  met  in  real  world  environments.  The  concept  of  discrete 
point  masses  is  not  very  realistic  either  as  far  as  geology  is 
concerned;  however,  thu  geodetic  goal  is  not  to  model  geological 
features  adequately  ^nd  accurately,  but  to  find  an  easy,  fast  and 
simple  way  to  model  the  external  anomalous  gravity  field;  point 
masses  present  themsiL-es  as  one  of  its  simplest  representations. 
Moreover,  using  point  masses  at  discrete  points  does  not  necessarily 
mean  that  masses  ar$-  concentrated  at  points,  since  the  gravity  field 
of  a  point  mass  is  not  distinguishable  from  that  of  a  spherically 
symmetric  mass  distribution.  Therefore,  we  could  as  well  think  in 
terms  of  a  continuous  anomalous  mass  distribution  which  is  such 
that  it  can  be  uniquely  described  by  a  point  mass  model.  In  this 
chapter  we  are  primarily  concerned  with  the  gravity  anomaly  co- 
variance  function  which  is  generated  by  a  discrete  distribution  of 
point  masses . 

Let  us  assume  that  I  point  masses  { y  ^ } ,  i  =  1,  ...,  I 
are  given  at  the  points  {Q^}  ,  i  =  1,  ...,  I  which  are  located  on 
a  concentric  sphere  with  radius  R-D  ;  we  assume  furthermore  that 
the  gravitational  constant  G  is  already  contained  in  u  .  The 
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Taking  into  account  the  representation  of  the  harmonic  function 
l~l  in  terms  of  a  Fourier  series  (Heiskanen  and  Moritz,  1967,  d.  33) 
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(3.4) 


with  a  =  1  -  ^  ,  and  using  the  orthonormality  relations  of  the 
fully  normalized  spherical  harmonics 
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(<5nr  ...  Kronecker-symbol) ,  the  integral  in  equation  (3.3)  1 

reduces  to 
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do (P)  = 
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and  the  Fourier  coefficients  assume  the  simple  form 


-  1  an  1  - 

anm  =  R  2n+l  . ^  U i * nm  (Qi 5 
i=l 

let  us  repeat:  anm  is  the  Fourier  coefficient  of  the  potential, 
generated  by  point  masses  { u i >  =  (y(Qi)  ,  i  =  1,  ...,  I 
Denoting  by 


W_L 

R 


(3.6a) 
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the  individual  potentials,  the  coefficients  a  are  given  by 


a 

run 


n 


2n+l  >  wi*nm(Qi) 

i  =  l 


(3.3)  " 


and  the  point  mass  generated  potential  by 


T  (P)  =  l 


n 


,  ,  “i  I  2 fer  ?  <v 

1=1  n=0  m  =  -ri 


(3.1)  ’ 


At  this  point  it  is  very  instructive  to  consider  the  special  case 
of  D  approaching  R  ;  in  other  words,  we  investigate  equation 
(3.1)'  for  the  extreme  case  of  all  point  masses  concentrated  at 
the  origin  of  the  coordinate  system.  As  a  matter  of  fact,  the 
corresponding  potential  should  reduce  to  that  of  a  single  centered 
mass  point 


GM 

R 


It  can  be  easily  shown  that 


lim  an  =  limj  1  -  ^ 

D-*-R  D-*-R  '• 


n ,  0 ' 


all  powers  of  a  vanish  apart  from  power  Zero;  therefore,  only 
ago  is  different  from  zero  and 
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T (P)  =  a00 


I 


GM 

R 


is  the  corresponding  potential  as  anticipated. 


Transition  to  gravity  anomalies. 

The  gravity  anomaly  field  at  sea  level  is  related  to  the 
anomalous  potential  through  a  convolution  (Stokes  and  inverse 
Stokes  formula) ;  the  corresponding  relation  in  the  frequency 
domain  is  given  by  the  well-known  product  (Heiskanen  and  Moritz, 
1967 ,  p.  97) 


-  nr  T„  •  ,3-7) 

Therefore,  the  Fourier  coefficients  b  of  the  gravity  anomaly 
Fourier  series 

-  ll  Enm7nm(p)  (3.8) 

run 

are  related  to  the  Fourier  coefficients  of  the  anomalous  potential 
through 


(3.9) 


and  with 
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u  u  . 

■i  Hi 


R  R‘ 


(3.6b) 


the  coefficients  are  obtained  by 


b  =  iri  )'  u  i  (Q.  )  . 

nm  2n+l  i  nm  i 

i=l 


(3.10) 


Consequently,  the  gravity  anomaly  at  r  =  R  is  represented  by 


I  _  ■=  n  __ 

ig(P)  =  )'  u,  l  — rr -  -  JB(P1  )  . 

’  -v  i  2n+l  -  nm  nm  i 

1=1  n=0  m=-n 


or,  considering  the  decomposition  formula  (Heiskanen  and  Moritz,  196 
p.  33) 


pri(c°SiJjp  )  =  JivTT  l  ^1P)frm(Q)  ' 


(3.11) 


m=-n 


even  simpler  by 


■ig(P)  ■  l  ^  l  ari(n-  l)Pn(cos^p  ) 
i=l  n=0  'vi 


(3.8)  ’ 
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The  gravity  anomaly  degree  variances. 

In  order  to  study  the  statistical  characteristics  of  tne 
gravity  anomaly  field,  generated  by  a  finite  number  of  point 
masses,  we  need  to  know  the  power  of  the  field  distributed  over 
all  frequencies  (degrees)  n  ;  for  a  specific  n  this  power  is 
denoted  degree  variance  of  degree  n  and  given  by  (Heiskanen  and 
Moritz,  1967,  p.  259) 
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r. 
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I 

m  =  -  n 


b2 

nm 


(3.12) 


With  b  given  by  (3.10),  the  degree  variances  can  be  written 
nm 

explicitly 
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or,  interchanging  the  sequence  of  summation. 


c  -  v1  i  - 

Cn  l  2n+lJ  ^ 


_  n 

I 

ra  =  -n 


i>  (Q.  )  i  (Q.) 
y  n  m  i  n  m  j 


The  last  sum  represents,  apart  from  the  factor  (2n+l)  ,  the  decom¬ 

position  of  the  Legendre  polynomial  Pn  into  fully  normalized 
spherical  harmonics  (eq.  (3.11));  therefore,  c  reduces  to 
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g2n (n-1)2 
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The  double  sum  equals  the  quadratic  form 
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with 

y  •  s  (p^/  ^2#  •••»  y.J  O.laa) 

and  the  symmetric  matrix 


A  : 
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I  Pn(cos*i 

i  . 
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1  )  ,  Pr,  (  COS1*'  1  2  )  ' 

**  t 


Pn(COSV; fI) 
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(3.15b) 


The  matrix  A  can  be  considered  a  covariance  matrix  derived 

n 

from  the  covariance  function  CU)  »  P  (cos^)  ,  with  elements 

n 

of  depending  on  the  mutual  spherical  distance  between  the 

individual  point  masses.  As  a  consequence,  An  has  only  posi¬ 
tive  and  zero  eigenvalues  for  n  >  0  and  as  a  matter  of  fact, 

am  — 

u  A  y  >0.  This  property  guarantees  the  non-negativity  of  the 
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{ }  which,  in  turn,  represents  a  homogeneous  and  isotropic  gravity 
anomaly  covariance  function,  derived  from  a  finite  set  of  point 
masses,  in  the  frequency  domain. 


c 

n 


a2n (n-1 ) 2  =T.  = 

2n+l -  u  V 


(3.13)  ' 


The  corresponding  covariance  function  of  gravity  anomalies,  at  zero 
level,  is  obtained  through  (Heiskanen  and  Moritz,  1967,  p.  256) 


w  /y  ^ 

C(K»)  -  l  uTAnuPn(cos-M  ;  (3.16) 

n  =  2 

its  spatial  extension  can  be  easily  derived  by  covariance  propa¬ 
gation  in  the  frequency  domain,  using  the  upward  continuation 
operator  applied  to  gravity  anomalies  (Heiskanen  and  Moritz,  1967, 
pp.  88-89)  , 


C(P,Q)  =  l 

n  =  2 


Q.2H  (n-1 )  2  (  r2  ln+2=T  = 

^1—  l^J  “  VPn(COBW 


(3.17) 


and  (3.17)  can  be  simplified  to  become 
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The  corresponding  spatial  covariance  function  of  the  anomalous 
potential  is  given  by 


K(P,Q)  = 
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(3.18) 


It  is  both  interesting  and  instructive  to  compare  the  gravity 
anomaly  degree  variances  (2.10)  with  (3.13)'  :  equation  (2.10)  has 
been  derived  from  a  white  noise  anomalous  mass  model  at  depth  D  , 
eq.  (3.13)'  from  an  (arbitrary)  discrete  distribution  of  a  finite 
number  of  point  masses  located  at  the  same  depth  D  .  The  common 
properties  are  a)  the  transition  from  mass  to  gravity,  and  b)  the 
upward  continuation  of  gravity;  the  combination  of  both  is  repre¬ 
sented  by  the  common  degree-dependent  factor  a2r*  (n-1 )  V  (2n+l )  ; 
the  most  striking  difference  is  the  degree-independent  multiplication 
factor 


f  4tG''2  Mo 


in  (2.10)  which  is  due  to  the  introduced  white  noise  model  of  the 
anomalous  mass  distribution,  and  the  degree-dependent  multiplication 
factor  ITTAnu  *-n  (3.13)’,  which  is  determined,  for  each  n  ,  by 
the  actual  point  masses.  It  can  be  shown  that  (3.13)'  converges 
to  (2.10)  if  the  number  of  point  masses  becomes  infinite  with  no 
correlation  between  neighboring  point  masses.  Consequently,  (2.10) 
is  a  special  limit  case  of  (3.13)'  and  this  is  why  we  investigated 
it  before.  (3.13)'  is  general  and  can  be  used  with  real  point  mass 
data. 
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4.  GRAVITY  ANOMALY  DEGREE  VARIANCES  DUE  TO  REGULARLY 
DISTRIBUTED  POINT  MASSES 

Gravity  anomaly  degree  variances  represent  the  power  of  the 
gravity  anomaly  field  broken  up  by  degree.  A  distribution  of  point 
masses  (at  a  certain  depth  D  )  yields  a  gravity  anomaly  field  at 
zero  level  which  can  be  statistically  described  by  its  degree 
variances.  For  low  to  moderately  high  degree  variances  are  known 
from  observations.  Any  physically  meaningful  point  mass  model  has 
therefore  to  meet  an  essential  requirement:  it  should  generate 
gravity  anomaly  degree  variances  which  are  close  to  the  "observed" 
degree  variance  model.  Equation  (3.13)'  relates  discrete  point 
masses  to  corresponding  gravity  anomaly  degree  variances. 

Given  a  discrete  point  mass  distribution,  the  very  problem 
of  calculating  degree  variances  consists  obviously  in  the  numerical 
calculations  of  the  quadratic  forms. 


This  harmless-looking  expression  poses  severe  problems  if  the 
point  mass  distribution  is  irregular  and  if  w  is  large:  the 
elements  of  each  matrix  An  are  Legendre  polynomials  of  degree 
n  ,  evaluated  for  the  mutual  spherical  distance  between  the  data 
(point  masses) , 

ai"’  *  pn(oOS*ij>  •• 


therefore,  the  calculation  of  each  element  a^  requires  one 
spherical  distance  calculation.  The  recurrence  relation  of 
Legendre  polynomials  (Heiskanen  and  Moritz,  1967,  p.23) 


pn(t)  *  tPn  -  “T  Pn  9(t) 

n  n  n-l  nn-2 


carries  over  to  the  matrices  A  , 

n  ' 


A  =  A  sA  -  —  A  , 

n  n  1  n- 1  n  n-  2  ' 


(4.1) 


where  "  sa  "  denotes  the  Hadamard  product  of  matrices , 


C  =  A  E  B  ,  =„-•„.  bi} 


Despite  of  the  pleasant  feature  of  An  (note  that  the  spherical 
distances  have  to  be  calculated  only  once  —  in  order  to  set  up 
the  matrix  A,^  )  ,  a  calculation  of  the  quadratic  form  uTAn“ 
becomes  prohibitive  for  irregularly  distributed  point  masses, 
even  if  its  number  I  is  as  low  as  a  few  hundred. 

Is  our  problem  surmountable?  Yes,  it  is.  The  picture  changes 
dramatically  if  we  assume  a  regular  distribution  of  the  point 
masses  at  the  grid  points  of  a  geographical  grid.  In  this  case  we 
can  take  advantage  of  the  structure  of  Ar  ,  transform  the  quadratic 
form  into  the  spectral  domain  using  fast  Fourier  transform  methods, 
and  apply  Parseval’s  theorem.  This  method  is  not  new  to  geodesy; 
the  interested  reader  will  find  a  detailed  treatment  of  such  kind 
of  problems  in  (Heller  et  al.,  1977)  and  (Colombo,  1979).  Therefore, 
we  restrict  ourselves  to  a  comprehensive  and  somewhat  simplified 
presentation  of  the  procedure. 
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Assume  J  point  masses  to  be  distributed  with  constant 
spacing  along  a  single  parallel.  The  corresponding  discrete 
Fourier  transform  matrix  F  is  given  by  (Heller  et  al.,  1977, 
p.  23) 


with  the  imaginary  unit  i  =  /-l  .  The  discrete  Fourier  trans¬ 
form  vector  of  the  vector  of  point  masses  will  be  denoted  by  X 
and  is  obtained  through  the  linear  transformation 

X  =  F*y  ,  (4.3) 

where  "  *  "  denotes  the  complex  conjugate  transpose.  As  a  matter 
of  fact  X  is  in  general  a  complex  vector?  it  consists  of  a  real 
(X°)  and  imaginary  (X1)  part 

X  =  X°  +  iX1 


with  elements 
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If  A  would  be  the  unit  matrix  I  ,  the  norm  u  u  can  easily 
n  * 

be  shown  to  be  equal  to  XX,  which  is  the  corresponding  norm 
in  the  frequency  domain. 
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y T  y  =  X*X  ;  (4.5) 

this  is  the  discrete  form  of  Parseval's  theorem  (Cooley  et  al., 
1967).  The  proof  is  straightforward:  we  introduce  (4.4)  in  (4.5) 
and  obtain 
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The  expression  between  the  parentheses  equals  6  ,  (  5  ... 
Kronecker  symbol)  according  to  the  orthogonality  relation  of 
exponentials  (Brigham,  1974,  p.99), 
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and  the  identity 
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follows  immediately. 
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Our  problem  is  a  bit  more  complicated:  the  matrix  is 

not  a  unit  matrix,  but  a  full  one,  whose  elements  depend  only  on 
the  spherical  distance  of  the  point  masses  (and  on  the  degree  n  ) . 
Due  to  the  geometry  of  the  data  pattern,  which  is  fully  represented 
in  A  ,  the  matrix  has  a  very  special  structure:  A  is  a  Toeplitz 
circulant  matrix;  (row  k  equals  row  k-1  shifted  one  element  to  the 
right,  with  the  last  element  of  row  k-l  wrapped  around  to  the 
first  place  in  row  k  .)  Therefore,  An  has  only  J  distinct 
elements;  (in  our  case  even  only  int[(J+l)/2]  due  to  the  de¬ 
pendence  of  the  spherical  distance.)  It  is  very  essential  for  the 
following  that  circulant  matrices  are  diagonalized  under  the  dis¬ 
crete  Fourier  transformation.  Denoting  the  discrete  Fourier  trans¬ 
form  of  A  by  the  diagonal  matrix  A  , 
n  n 

\n  =  FA  F*  ,  (4.7) 

n  n 

the  diagonal  terms  are  simply  given  by 
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denote 

the  elements  of 

the  first 

row  of  An  .  The  proof  of  equation  (4.8)  is  simple;  it  relies  on 
the  circulant  property  of  A^  and  on  the  orthogonality  relation 
(4.6).  The  discrete  form  of  Parseval's  theorem  with  An  as 
metric  is  given  by 

X* A  X  =  mTA  Z  • 
n  n 


(4.9) 
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(Note  that  is  the  diagonal  matrix  of  eigenvalues,  F  the 

orthogonal  matrix  of  eigenvectors  of  A^  ;  X  is  the  image  of 
u  in  the  eigenvector  system.)  The  discrete  Fourier  transformations 
(4.3)  and  (4.7)  are  accomplished  by  the  enormously  powerful  tool 
of  Fast  Fourier  Transform  (FFT)  (Brigham,  1974;  Singleton,  R.C., 
1968).  (The  CPU-time  increases  with  J • InJ  ;  the  constant  multi¬ 
plication  factor  is  about  2.7  •  10  5  sec  for  the  Singleton-algorithm 
on  a  UNIVAC  1100/81.)  In  view  of  this  speed,  the  transformation 
into  the  frequency  domain  does  not  pose  any  problems.  The  matrix 
An  can  be  computed  recursively  using  the  relation  (4.1);  (actually 
only  half  of  the  elements  of  the  first  row  have  to  be  calculated 

because  of  the  circulant  character  of  A  .)  Parseval ' s  theorem 

n 

(4.9)  gives  immediately  the  desired  degree  variances  (apart  from 
the  factor  a 2n (n-1 ) 2 ( 2n+l ) ~ 1  ). 

Sofar  we  have  considered  a  regular  data  distribution  restricted 
to  a  single  parallel.  Let  us  now  investigate  the  case  of  a  regular 
distribution  on  K  parallels.  Now  the  data  vector  Z  consists  of 

5S  (  V  ) 

K  subvectors  y  with  J  elements  each.  Analogously,  the 

2 

symmetric  matrix  A  can  be  subdivided  into  K  submatrices  of 

2  °  f  v  3c '  ) 

dimension  J  each.  Each  submatrix  An  '  ,  which  corresponds 

to  row  k  combined  with  k’  ,  is  Toeplitz  circulant  due  to  the 
geometry  of  the  data  pattern.  Therefore,  A^  is  a  symmetric  block 
matrix  with  Toeplitz  circulant  blocks.  If  the  data  pattern  is 

symmetric  with  respect  to  the  equator,  A  is  even  persymmetric . 

=  ( k )  n  ( k  k  ’  ) 

Transforming  each  subvector  y  and  each  submatrix  A ' 

into  the  frequency  domain  according  to  equ.(4.4)  and  (4.8)  leads 

to  Parseval’ s  identity  with  a  block-diagonal  eigenvalue  matrix 

A  ,  corresponding  to  the  block-circulant  matrix  A  .  As  a  matter 
n  n 

of  fact,  a  persymmetry  of  An  is  also  reflected  by  A  .  More  de¬ 
tails  on  that  particular  procedure  are  contained  in  (Colombo,  1979) . 
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5.  MULTI-LAYER  MASS  DISTRIBUTIONS 


We  have  seen  in  chapter  2  that  a  single- layer  anomalous 
mass  distribution  is  neither  able  of  modelling  the  observed 
Earth's  gravity  field  adequately,  nor  is  it  very  realistic  from 
a  geophysical  point  of  view.  Therefore,  multi-layer  models  have 
been  proposed  by  numerous  authors  (Schwarz,  1977,  1981;  Jordan, 
1978;  Heikkinen,  1981;  Davenport,  1978)  .  In  the  sequel  we  in¬ 
vestigate  the  multi-layer  mass  distribution  problem  with  an 
emphasis  on  the  statistics  of  the  resulting  anomalous  gravity 
field.  In  the  same  order  as  before,  we  consider  first  white 
noise  distributions  and  second  discrete  point  mass  models  with 
regular  distribution. 

Denoting  the  depths  of  the  mass  anomaly  layers  by  D 
(  1  =  1 , . . . ,L  ),  the  eigenvalues  of  the  integral  kernel,  which 
is  responsible  for  the  transition  from  mass  to  gravity,  are  ob¬ 
tained  from  (2.7)  , 
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(5.1) 


with  a.  =  1  -  .  The  gravity  anomaly  power  spectrum  corre- 

1  K 

sponding  to  the  white  noise  anomalous  mass  distribution  at  depth 
Dl  is  obtained  by  (2.10) 
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4itG 


2  M 


(1) 
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(n-1)2  2n 
2n+l  °i 


(5.2) 


and  the  total  power  spectrum  is  obviously  the  sum  of  the  indi¬ 
vidual  ones. 
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provided  two  mutually  different  mass  layers  are  uncorrelated.  The 
total  gravity  anomaly  covariance  function  is  formally  equal  to 
(2.11)  with  gn  defined  by  (5.3), 


C 

gg 


(t) 


l  ?nPn(t) 
n=2  n  n 


If  the  anomalous  mass  distribution  is  discrete  and  regular, 
such  that  the  data  pattern  is  common  to  all  levels,  the  disturbing 
potential  harmonic  coefficients  (3.3)"  are  represented  in  terms  of 


1 

2n+l 


1  _  L 
y  *  (q  )  y  u  an  , 

.Vn“  i  ,  ,  il  1 
1=1 


where  u  denotes  the  point  mass  (multiplied  by  the  gravitational 
constant  G  and  divided  by  the  mean  Earth  radius  R  ) ,  which  has 
Q.^  as  horizontal  and  R-D^  as  vertical  position.  (Note  that 
can  be  considered  as  the  potential,  which  is  generated  by 
that  particular  point  mass  on  a  ball  with  radius  R  and  centered 
at  the  point  mass.)  Analogously,  the  gravity  anomaly  coefficients 
(3.10)  are  given  by 


b 

run 


2n+l  .  E  pnm(Qi)  J-  ^  i  1 a  1  * 
i=l  1=1 


(5.5) 


For  a  specific  degree  n  the  last  sum  represents  (apart  from  the 
2 

constant  G/R  )  the  weighted  sum  of  mass  points,  which  are  located 
on  the  same  radius  vector;  the  weights  are  given  by  the  individual 
depths.  Denoting  this  sum  by  , 
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( n ) 
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I 


1  = 


n 

I 


(  5 . 6  ' 


the  gravity  anomaly  degree  variances  c  can  be  represented  in 
a  similar  form  as  in  equation  (3.13)', 


c 


n 


( n- 1 )  ^  =  (n)  T,  =  ( n ) 

2nTl  '  u  V 


(5.7) 


All  degree  variances  (for  arbitrarily  high  n  )  are  theoretically 
affected  by  all  point  masses.  However,  due  to  the  weight  factors 
a"  ,  deep  point  masses  have  practically  no  significant  contribution 
to  high  degree  variances.  Needless  to  say,  the  gravity  anomaly 
variance  equals  the  sum  of  all  degree  variances. 
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